CP2 skyrmions and skyrmion crystals in realistic quantum magnets

Magnetic skyrmions are nanoscale topological textures that have been recently observed in different families of quantum magnets. These objects are called CP1 skyrmions because they are built from dipoles—the target manifold is the 1D complex projective space, CP1 ≅ S2. Here we report the emergence of magnetic CP2 skyrmions in a realistic spin-1 model, which includes both dipole and quadrupole moments. Unlike CP1 skyrmions, CP2 skyrmions can also arise as metastable textures of quantum paramagnets, opening a new road to discover emergent topological solitons in non-magnetic materials. The quantum phase diagram of the spin-1 model also includes magnetic field-induced CP2 skyrmion crystals that can be detected with regular momentum- (diffraction) and real-space (Lorentz transmission electron microscopy) experimental techniques.

Lord Kelvin's vision of the atom as a vortex in ether 1 inspired Skyrme 2,3 to explain the origin of nucleons as emergent topologically non-trivial configurations of a pion field described by a 3 + 1 dimensional O(4) non-linear σ-model. In the modern language, these "skyrmions" are examples of topological solitons, and Skyrme's model has become the prototype of a classical theory that supports these solutions. Besides its important role in high-energy physics and cosmology, Skyrme's model also led to important developments in other areas of physics. For instance, the baby Skyrme model 4-6 (planar reduction of the nonlinear σ-model), which is an extension of the Heisenberg model 4,5,7 , has baby skyrmion solutions in the presence of a chiral symmetry-breaking Dzyaloshinskii-Moriya (DM) interaction [8][9][10][11] .
The target manifold of the above-mentioned planar baby skyrmions is S 2 ≅ CP 1 , i.e., the usual 2D sphere, corresponding to normalized dipoles. More generally, one may consider the complex projective space CP N−1 that represents the normalized N-component complex vectors, up to an irrelevant complex phase. The topologically distinct, smooth mappings from the base manifold S 2 (2D sphere ≅ compactified plane) to the target manifold CP N−1 can be labeled by the integers: Π 2 ðCP NÀ1 Þ = Z. This homotopy group suggests generalizations of the planar Skyrme's model to N > 2, such as the CP 2 nonlinear σ-model [31][32][33] and in the Faddeev-Skyrme type model 34,35 . In condensed matter physics, chiral CP 2 skyrmion configurations induced by fluctuations or quenching the system through a phase transition were proposed in the context of three-band superconductors with broken time-reversal symmetry [36][37][38] . In recent work, Akagi et al. considered the SU(3) version of the Heisenberg model with a DM interaction, whose continuum limit becomes a gauged CP 2 nonlinear σ-model with a background uniform gauge field 39 . An attractive aspect of this model is that it admits analytical solutions by the application of techniques developed for the gauged non-linear σ-model. However, it may be challenging to find materials described by this model because SU (3) can only be an accidental symmetry of the spin-spin interactions of real quantum magnets, and Hamiltonians that do exhibit SU(3)-invariance contain unrealistically strong biquadratic terms. In insulating magnets, biquadratic interactions are typically much smaller than bilinear interactions because they are of higher order in the small parameter that leads to the emergence of magnetic moments (localized electrons) in real materials (e.g., the ratio t/U between the typical hopping amplitude, t, and the on-site Hubbard repulsion, U, in the case of Mott insulators). Similar limitations apply to other works that study skyrmion solutions of the bilinear-biquadratic spin one model [40][41][42][43] .
The main purpose of this work is to demonstrate that exotic CP 2 skyrmions readily emerge in a simple and realistic spin-1 (N = 3) model and its natural extensions. In other words, we propose that these magnetic textures could likely be observed in real materials. Remarkably, isolated CP 2 skyrmions can either be metastable states of a quantum paramagnet (QPM) or a fully polarized (FP) ferromagnet. Unlike the "usual" CP 1 magnetic skyrmions, the dipolar field of the metastable CP 2 skyrmions of quantum paramagnets vanishes away from the skyrmion core. Moreover, the application of an external magnetic field to the QPM induces stable triangular crystals of CP 2 skyrmions in the field interval that separates the QPM from the FP state.

Model
To illustrate the basic ideas, we consider a minimal spin-1 model defined on the triangular lattice (TL): The first term includes an easy-axis ferromagnetic (FM) nearestneighbor exchange interaction J 1 < 0 and a second-nearest-neighbor antiferromagnetic (AFM) exchange J 2 > 0. For simplicity, we assume that the exchange anisotropy, defined by the dimensionless parameter Δ > 1, is the same for both interactions. The second and third terms represent the Zeeman coupling to an external field and an easy-plane single-ion anisotropy (D > 0).Ĥ is invariant under the space group of the TL and the U(1) group of global spin rotations along the field axis. We will adopt |J 1 | as the unit of energy (i.e. J 1 = −1).
To study the properties of the skyrmion solutions of Eq. (1), it is helpful to consider the classical limit first. In doing so, we follow the existing literature on topological solitons, which are inherently classical objects. There is, moreover, good reason to expect that quantum fluctuations are not relevant to the present study. Experimentally, there is existing evidence of spin-1 triangular materials that exhibit semi-classical spiral orderings due to competing ferromagnetic and antiferromagnetic exchange interactions 44 . Some of these are discussed in Section "Disussion". Furthermore, the results that will be developed will remain unaltered for a simple 3D extension of the current model, achieved by vertically stacking triangular layers with ferromagnetic interlayer coupling. The larger coordination number of the 3D model and the long wavelength nature of the ordered states both act to reduce quantum fluctuations, further justifying the classical approximation.
It is important to note, however, that there are subtleties in formulating the appropriate classical limit 45,46 . The traditional classical limit is based on SU(2) coherent states, which retain only the spin dipole expectation value and produces the Landau-Lifshitz spin dynamics. This approach is adequate for modeling systems with weak single-ion anisotropy D ≪ |J 1 |. To classically model systems in the regime D ≳ |J 1 |, however, it is necessary to retain more structure from the quantum spin-1 states, which live in a local Hilbert space of dimension N = 3. Specifically, our classical limit will assume that the many-body quantum state is a direct product of SU(3) coherent states [45][46][47][48][49][50][51][52] : where Z j = ðZ 1 j ,Z 2 j ,Z 3 j Þ T is a complex vector of unit length and f|x 1 i j ,|x 2 i j ,|x 3 i j g is an orthonormal basis for the local Hilbert state onsite j.
Local physical operators are represented by Hermitian matrices that act on SU(3) coherent states. The space of 3 × 3 traceless, Hermitian matrices comprises the fundamental representation of the suð3Þ Lie algebra. A basisT where we are using the convention of summation over repeated Greek indices. We may additionally impose an orthonormality condition It is conventional to define the structure constants as acting on site j are generators for a spin-1 representation of SU (2). It is possible to formulate generators of SU(3) as polynomials of these spin operators, where T 7,5,2 j are the dipolar components of the spin-1 degree of freedom, while the other five generators are the quadrupolar components. Here we have adopted the notation and conventions of ref. 39 to make closer contact with the literature on high-energy physics. (Our definitions forŜ x andŜ z differ from these two in ref. 39 by a minus sign). Let |1i j , |0i j , and | 1 j denote the normalized eigenstates ofŜ z j , with eigenvalues, 1, 0 and −1, respectively. In the Cartesian basis, the SU(3) generators given in Eq. (5) are the Gell-Mann matrices: The orbit of coherent states |Z j i is obtained by applying SU(3) transformations to the highest weight state |1i j 45 : |Z j i =Û j |1i j . Since the global phase is a gauge degree of freedom, the orbit of physical SU(3) coherent states is S 5 /S 1 ≅ CP 2 . The "SU(3) classical limit" of the spin Hamiltonian (1) is obtained by replacing the Hamiltonian operatorĤ with its expectation value after rewritingĤ in terms of SU(3) spin components, where I μ ij = J ij ðΔδ μ,2 + δ μ,5 + δ μ,7 Þ and B μ = ðÀhδ μ,2 À Dδ μ,8 = ffiffiffi 3 p Þ. Because the direct product form of Eq. (2), H can be expressed as a function of the "color field" which satisfies the constraints where d μνη = 1 4 Tr ðλ μ fλ ν ,λ η gÞ. This in turn leads to the Casimir identity: d mpq n m n p n q = 8 9 : In terms of this color field, we can express To avoid an explicit use of the structure constants ( f ημν ), we introduce an equivalent formulation of the problem using the operator field n j = n μ j λ μ . Topological soliton solutions of the color field become well-defined in the continuum (long wavelength) limit, where the Hamiltonian can be approximated by where ∇ denotes the spatial gradient operator. The coupling constants can be expressed in terms of the parameters of the lattice model (9): Eq. (13) corresponds to an anisotropic CP 2 model. For skyrmion solutions, the base plane R 2 can be compactified to S 2 because the color field takes a constant value n ∞ at spatial infinity. These spin textures can then be characterized by the topological charge of the mapping n : R 2 ∼ S 2 7 !CP 2 : For the lattice systems of interest, the CP 2 skyrmion charge can be computed after interpolating the color fields on nearest-neighbor sites n j and n k along the CP 2 geodesic: where △ jkl denotes each oriented triangular plaquette of nearestneighbor sites j → k → l and γ kj = arg½hZ k |Z j i is the Berry connection on the bond j → k (see Supplemental Material). We emphasize that the color field formalism just discussed is fully equivalent to the formalism based on coherent states. In particular, it is straightforward to show that the operator representation of the color field may be expressed as n j = |Z j ihZ j | À 1=3, in which form it becomes clear that both n and the coherent state |Z j i provide equivalent representations of the same classical state 53,54 .

Phase diagram
The T = 0 phase diagram (Fig. 1) is obtained by numerically minimizing the classical spin Hamiltonian H (12) in the 4L 2 -dimensional phase space of a magnetic cell of L × L spins (see the "Methods" section). The shape and size of this unit cell is dictated by the symmetry-related magnetic ordering wave vectors Q ν (ν = 1, 2, 3) (see Fig. 2a, b), which are determined by minimizing the exchange interaction in momentum space: JðqÞ = P jl J jl e iqÁðr j Àr l Þ : The ratio between both exchange interactions, J 2 =| J 1 | = 2=ð1 + ffiffi ffi 5 p Þ is tuned to fix the magnitude of the ordering wave vectors, |Q ν | = |b 1 |/5 27 , corresponding to a magnetic unit cell of linear size L = 5. As we will see later, the relevant qualitative aspects of the phase diagram do not depend on the particular choice of the model (see the section "Large-D limit"). The three ordering wave vectors, which are related by the C 6 symmetry of the TL, are parallel to the Γ-M ν directions (denoted in Fig. 2).
The resulting phase diagram shown in Fig. 1 includes multiple magnetically ordered phases between the FP phase and the QPM phase, where every spin is in the |0i state. For D ≫ |J 1 |, these phases include two field-induced CP 2 skyrmion crystals (SkX-I and SkX-II), separated by two modulated vertical spiral phases (MVS-I and MVS-II), where the polarization plane of the spiral is parallel to the c-axis and the magnitude of the dipole moment is continuously suppressed as the moment rotates fromẑ to Àẑ directions. The spiral phases have the same symmetry and are separated by a first-order metamagnetic transition. As shown in Fig. 2a, the CP 2 skyrmions of the SkX-I crystal have dipole moments that evolve continuously into the purely nematic state (|0i) as they move away from the core. Conversely, Fig. 2b shows that the spins in the SkX-II phase have a strong quadrupolar character (the small dipolar moment is completely suppressed in the large D/|J 1 | limit) at the skyrmion core, and evolve continuously into the magnetic state |1i as they move away from the core. The CP 2 skyrmion density distribution ρ jkl is also indicated with colored triangular plaquettes in Fig. 2a, b for SkX-I and SkX-II, respectively. As shown in the inset of Fig. 1, phase SkX-II extends down to D/| J 1 | ≃ 5, while phase SkX-I disappears near D/| J 1 | ≃ 8.
New competing orderings appear in the intermediate D/|J 1 | region. In particular, a significant fraction of the phase diagram is occupied by the so-called canted spiral (CS) phase, |Z j i = cos θ|0i j + e iQÁr j sin θ cos ϕ|1i j + e ÀiQÁr j sin θ sin ϕ| 1 j , ð17Þ where θ and ϕ are variational parameters, and Q can take any values among {Q 1 , Q 2 , Q 3 }. Upon increasing D, the magnitude of the dipole moment of each spin, |hŜ j i|, is continuously suppressed to zero at the boundary, that signals the second-order transition into the QPM phase. As shown in Fig. 1, several competing phases appear above the CS phase upon increasing h. These phases include three triple-Q spiral orderings [3Q I-III] with dominant weight in one of three Q transverse components and a staggered distribution of the CP 2 skyrmion density ρ jkl [see Eq. (16)] and three different modulated double-Q orderings (MDQ I-III) and two triple-Q orderings. All of these phases are described in detail in the supplementary information. In the rest of the paper, we will focus on the SkX phases and the single-skyrmion metastable solutions that emerge in their proximity.

Large-D limit
The origin of the CP 2 skyrmion crystals can be understood by analyzing the small |J ij |/D regime, whereĤ can be reduced via first-order degenerate perturbation theory in J ij /D to an effective pseudo-spin-1/2 low-energy Hamiltonian, The pseudo-spin-1/2 operators are the projection of the original spin operators into the low-energy subspace S 0 generated by the quasidegenerate doublet f|0i j ,|1i j g (see Fig. 3): where P 0 is the projection operator of the low-energy subspace. Importantly, the first state of the doublet has a net quadrupolar moment but no net dipole moment, 0 h |Ŝ j |0i j = 0, while the second state maximizes the dipole moment along theẑ-direction 1 h |Ŝ j |1i j =ẑ. This means that three pseudo-spin operators generate an SU(2) subgroup of SU(3) different from the SU(2) subgroup of spin rotations.
H eff represents an effective triangular easy-axis XXZ model with effective exchange, anisotropy and field parametersJ ij = 2J ij ,Δ = Δ 2 and h = h À D À 3ΔðJ 1 + J 2 Þ, respectively. This model is known to exhibit a field-induced CP 1 SkX phase 25, 27 on a lattice for fixed choice of J 1 and J 2 . Further study has demonstrated that the full field-anisotropy phase diagram remains qualitatively the same upon approaching the long wavelength limit of J 2 ! 1 3 |J 1 |, the Lifshitz point where the ordering wave vectors go to zero 26 . It follows that lattice effects do not alter the qualitative features of the phase diagram for wavelengths at least as long as that set by the J 1 and J 2 examined here. In other words, these results do not depend on a fine-tuning of exchange parameters. Indeed, the continuum model for Eq. (20) matches the universal Hamiltonian presented in 26 , where η = x, y, z denotes the three components of the unit vector field n (|ñ| = 1), and where s = 1/2. Although the target manifold of this theory is CP 1 (orbit of SU(2) coherent states that belong S 0 ), we must keep in mind that H eff describes the large D/|J 1 | limit where the CP 2 skyrmions of the . The color scale of the arrows indicates hŜ z j i = À n 2 j . The insets display the static spin structure factors S ? ðqÞ = n 7 q n 7 Àq + n 5 q n 5 Àq and S zz ðqÞ = n 2 q n 2 Àq , with n q = P j e iqÁr j n j =L. The CP 2 skyrmion density distribution ρ jkl [see Eq. (16)] is indicated by the color of the triangular plaquettes. original spin-1 model become asymptotically close to CP 1 pseudo-spin skyrmions. In other words, the SkXs include a finite | 1 component for finite D/|J 1 |, which increases upon decreasing D. This component, which only appears in the low-energy model when second-order corrections in J ij are included, is responsible for the metamagnetic transition between the MVS-I and MVS-II phases (the transition disappears in the D → ∞ limit).
SinceĤ eff ðhÞ andĤ eff ðÀhÞ are related by a pseudo-time-reversal (PTR) transformation (ŝ j ! Àŝ j on the lattice andñ ! Àñ in the continuum) their corresponding ground states are related by the same transformation. In particular, the ground state ðñ =ẑÞ that is obtained above the saturation field (B >B sat ) corresponds to the FP state (hŜ j i =ẑ) in the original spin-1 variables, while the ground state ðñ = À zÞ below the negative saturation field (B< ÀB sat ) corresponds to the QPM phase (|Z j i = |0i j ). Correspondingly, the SkX induced by a positive h has pseudo-spins polarized along the quadrupolar direction (|0i) near the core of the skyrmions and parallel to the dipolar one (|1i) at the midpoints between two cores. This explains the origin of the SkX-II crystals depicted in Fig. 2b. The negative B counterpart of this phase, which is obtained by applying the PTR transformation, corresponds to the SkX-I crystal shown in Fig. 2a. In this case the skyrmion cores are magnetic, while the midpoints are practically quadrupolar (they become purely quadrupolar in the large D/|J 1 | limit). This simple reasoning explains the origin of the novel SkX phases included in the T = 0 phase diagram of H shown in Fig. 1. The intermediate phase between the SkX-I and SkX-II ground state of H induced by positive and negative values of h is a single-Q spiral with a polarization plane parallel to the caxis known as a vertical spiral (VS). This explains the origin of the MVS-I and MVS-II phases in between the two SkX phases (the first-order transition between both phases disappears in the large-D limit 25 ).

Single-skyrmion solutions
Besides the SkX phases shown in Fig. 4, the effective field theory (21) is known to support metastable CP 1 single-skyrmion solutions beyond the saturation fields |B| >B sat . The pseudo-spin variable is anti-parallel to the external field at the core and it gradually rotates towards the direction parallel to the field upon moving away from the center. Interestingly, this pseudo-spin texture translates into metastable single-skyrmion solutions of the QPM phase that have a magnetic core and a nematic periphery, as it is shown in Fig. 4a and b for different sets of Hamiltonian parameters. The CP 2 skyrmions are metastable solutions in the QPM phase for D ≳ 14, implying that these exotic magnetic-nematic textures should emerge in real magnets under quite general conditions. Similarly, the metastable pseudo-spin single-skyrmion solutions of the FP phase (B >B sat ) lead to a spin texture with a nematic (nonmagnetic) core and a magnetic (FP) periphery, like the one shown in Fig. 4c. Interestingly, this exotic CP 2 skyrmion solution remains metastable down to D ≃ 4|J 1 | and it coexists with regular (CP 1 ) metastable skyrmion solutions, like the one shown Fig. 4d, that emerge below D ≃ 4.25|J 1 |.

Discussion
We have demonstrated that CP 2 skyrmion textures emerge in realistic models of hexagonal magnets out of the combination of competing exchange interactions and single-ion anisotropy. It is important to note that the skyrmion crystals and metastable solutions reported in this work survive in the long wavelength limit 26 , implying that the CP 2 skyrmion phases described here should also exist in extensions of the model to honeycomb and Kagome lattice geometries.
There are a number of candidate materials that are well described by the spin-1 model given in Eq. (1). In particular, one may point to the series of triangular antiferromagnets of the form of ABX 3 , BX 2 , and ABO 2 44,55,56 , where A is an alkali metal, B is a transition metal, and X is a halogen atom. Compounds, such as FeI 2 57,58 , are described by the Hamiltonian of Eq. (1), but the sign of the single-ion and exchange anisotropies is opposite to the case of interest in this work. Related compounds, such as CsFeCl 3 , are known to be quantum paramagnets described by the same Hamiltonian with a dominant easy-plane singleion anisotropy D/J 1 ≃ 10 59 . An alternative route to finding realizations of our spin-1 Hamiltonian is to consider hexagonal materials comprising 4f magnetic ions with a singlet single-ion ground state and an excited Ising-like doublet (the effective easy-plane single-ion anisotropy D is equal to the singlet-doublet gap). Ultracold atoms are also powerful platforms to realize spin-1 models with tunable single-ion anisotropy 60 .
While a full examination of the new response functions and functionalities of the CP 2 skyrmions must be left to future research, a few remarks should be made here. It is clear that the intrinsically inhomogeneous nature of the local order parameter, which evolves from dipolar to quadrupolar upon moving toward or away from the skyrmion core, can lead to new behaviors. For instance, metastable CP 2 skyrmions above the saturation field can become stable (ground state) solutions by increasing the D term of a given magnetic ion. This can be achieved with the insertion of non-magnetic impurities that modify the local crystal field. Correspondingly, it should be possible to induce metastable CP 2 skyrmions by dynamically varying the local crystal field that determines the value of D. Furthermore, CP 2 skyrmions can be manipulated by applying a local strain due to the characteristically non-uniform distribution of the magnitude of their quadrupolar moment.
Before concluding, we remark on a subtle mathematical point. By definition, CP 2 skyrmions are distinguished from the more familiar CP 1 skyrmions by their enlarged target manifold. This distinction can be physically relevant: a CP 2 skyrmion will typically have a combination of dipolar and quadrupolar structures. The presence of quadrupole degrees of freedom will bring additional dynamical modes and will have entropic consequences. From a topological perspective, however, there is a certain sense in which CP 1 and CP 2 skyrmions are equivalent. To elaborate on this point, we first remark that CP 1 is a submanifold of CP 2 . Further, any CP 1 skyrmion can be faithfully embedded in the space of CP 2 skyrmions, and this embedding preserves skyrmion winding number. Such an embedded spin texture can be smoothly deformed to any other CP 2 skyrmion with an equal winding number, which establishes a topological equivalence.
In summary, this paper demonstrates that novel magnetic fieldinduced CP 2 skyrmion crystals should emerge in the presence of competing ferromagnetic and antiferromagnetic exchange interactions, a moderate easy-axis exchange anisotropy Δ > 2, and a dominant single-ion easy-plane anisotropy D that is strong enough to stabilize a QPM at T = 0. The field-induced quantum phase transition between the uniform quadrupolar state induced by the strong single-ion anisotropy and the CP 2 skyrmion crystal is presaged by the emergence of metastable CP 2 single-skyrmion solutions exhibiting a magnetic skyrmion core that decays continuously into a quadrupolar periphery. These novel skyrmions can be induced by applying a sufficiently large magnetic field to quantum paramagnets with competing exchange interactions and they can be manipulated with local strain.
The general principles discussed in this work can be generalized to N-level systems to obtain CP N−1 skyrmion crystals solutions from realistic spin Hamiltonians, illustrating the rich diversity of topological textures that can emerge in magnetic materials due to the quantum mechanical nature of their magnetic moments.

Methods
The numerical minimization for the phase diagram Fig. 1 is done in a cell of 10 × 10 spins containing four magnetic unit cells (L = 5). Two crucial steps are useful to improve the efficiency of the local gradientbased minimization algorithms 61 . In the first step, we set multiple random initial conditions |Zi (~50 for our case), where |Z j i on every site j is uniformly sampled on the CP 2 ≃ S 5 /S 1 manifold. After running the minimization algorithm, we keep the solution with the lowest energy for a given set of Hamiltonian parameters. In the next step, half of the initial conditions are randomly generated, while the other half corresponds to the lowest-energy solutions obtained in the first step within a predefined neighborhood of the Hamiltonian parameters. This procedure is iterated until the phase diagram converges.

Data availability
All data presented in this study can be reproduced using the code package described in the section "Code availability".

Code availability
The algorithms used in our numerical simulations are described in the "Methods" section. The numerical code is implemented in Julia and can be found at https://github.com/Hao-Phys/CP2Skyrmions.jl.